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The space based gravitational wave detector LISA is expected to observe a large population of Galactic white 
dwarf binaries whose collective signal is likely to dominate instrumental noise at observational frequencies in 
the range IIP 4 to IIP 3 Hz. The motion of LISA modulates the signal of each binary in both frequency and 
amplitude, the exact modulation depending on the source direction and frequency. Starting with the observed 
response of one LISA interferometer and assuming only doppler modulation due to the orbital motion of LISA, 
we show how the distribution of the entire binary population in frequency and sky position can be reconstructed 
using a tomographic approach. The method is linear and the reconstruction of a delta function distribution, 
corresponding to an isolated binary, yields a point spread function (psf). An arbitrary distribution and its recon¬ 
struction are related via smoothing with this psf. Exploratory results are reported demonstrating the recovery of 
binary sources, in the presence of white Gaussian noise. 


I. INTRODUCTION 

The Laser Interferometric Space Antenna (LISA) jlll is a 
space based interferometric gravitational wave (GW) detec¬ 
tor, currently in design, with a launch date sometime around 
2015. LISA will consist of three proof masses in drag free or¬ 
bits around the Sun. Passing GW signals will be detected by 
monitoring fluctuations in the nominal distance of 5 x 10 6 km 
between the proof masses using laser interferometry. LISA 
will be most sensitive to GW signals in the 10 4 to 1 Hz band 
and is expected to observe a wide variety of GW sources dur¬ 
ing its 3 to 5 year mission liftime. 

Among the sources predicted for LISA will be close white 
dwarf binaries within the Galaxy, estimated to number around 
10 8 or more in the low frequency range (10 4 to 10 3 Hz) of 
the LISA sensitivity band y]. LISA will be an all-sky moni¬ 
tor and the signals from long-lived GW sources from the en¬ 
tire sky will be summed in the response of the detector. The 
sum of the signals from all the Galactic binaries, which we 
call the binary population signal, is expected to dominate the 
instrumental noise of LISA at low frequencies and may act 
as a noise source from which signals from other sources, in¬ 
cluding bright individual binaries, will need to be extracted. 
On the other hand, the spatial and period distribution of LISA 
binaries is in itself an important astrophysical observable for 
the LISA mission. Measuring the distribution function would 
yield information pertinent to studies of the Galactic structure, 
for instance. 

The reconstruction of the distribution function of sources, 
given the superposed signal from all of them, is a critical data 
analysis problem for LISA. Intertwined with this task is that 
of removing as much of the binary population signal as pos¬ 
sible from the data so that other types of sources in the same 
frequency range can be detected. Note that the binary distri¬ 
bution function automatically includes the case of an isolated 
binary, in which case it is a delta-function. Hence, reconstruc¬ 
tion of the distribution function can also identify well isolated 
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individual binaries. Several methods Hill are currently 
being explored in the above context. 

In this paper, we present the first steps of a promising ap¬ 
proach to reconstructing the LISA binary distribution func¬ 
tion. The approach, which we call tomographic reconstruc¬ 
tion , is based on the observation that, loosely speaking, the 
response of LISA to the binary population signal is a Radon 
transform @] of the true binary distribution function. The ori¬ 
gin of this lies in the doppler modulation of binary signals due 
to the orbital motion of LISA. Given the measured response 
of LISA it is, therefore, possible to reconstruct the distribu¬ 
tion function by applying an inverse Radon transform to the 
response. 

The main purpose of the paper is to demonstrate the con¬ 
nection between the response of LISA and the Radon trans¬ 
form of the binary distribution function. We then explore 
how one of the most straightforward inversion methods for 
the Radon transform, namely direct Fourier inversion, can be 
used to reconstruct the binary distribution function. Tomo¬ 
graphic reconstruction is a linear operation on the data. Thus, 
it is fully characterized by the response to a delta-function dis¬ 
tribution or, in other words, by the point spread function (psf). 
The reconstructed distribution is related to the true one via 
smoothing by this psf. 

We present exploratory results using a simple, and cer¬ 
tainly naive, method to identify source candidates in the re¬ 
constructed distribution function. We also use a much simpli¬ 
fied scenario for demonstrating reconstruction. For instance, 
we only take the doppler modulation due to the orbital motion 
of LISA into account. We find that, for an observation pe¬ 
riod of 1 year and frequencies near 2 mHz, binaries having a 
matched filtering signal to noise ratio (SNR) of ~ 7 or higher 
and a minimum spacing of about 10 frequency bins (1 bin 
~ 3 x 10“ 8 Hz) can be distinguished fairly easily. This is in 
the range of the expected density of bright binaries d. Bina¬ 
ries can also be seen at separations of 4 to 6 frequency bins 
but the simple identification method finds spurious sources 
arising from the superposition of the sidelobes of the psf. A 
more sophisticated identification method based on the use of 
deconvolution methods, such as CLEAN d , is required in 
this case. However, we will address the latter issue in a future 
work since it is a fairly non-trivial problem in itself. 

The rest of the paper is organized as follows. In Section Hfl 
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we give an outline of the Radon transform and the Direct 
Fourier Inversion method. Section mu describes how the re¬ 
sponse to a distribution function of binaries can be expressed 
in terms of a Radon transform. The algorithm for tomographic 
reconstruction is presented in Section ITvl SectionlVlcontains 
exploratory results. 


H. RADON TRANSFORM AND ITS INVERSION 

The Radon transform @] of a function /(x), x € R”, is de¬ 
fined in terms of the integrals of /(x) over all possible hyper¬ 
planes of dimension n — 1. The function /(x) can be recon¬ 
structed completely from its Radon transform provided inte¬ 
grals are available for all hyperplanes in the original domain. 

An important application area of Radon transform is in 
medical imaging where it is well known as the CAT scan. 
The absorption of X-rays passing through the three dimen¬ 
sional density distribution of biological matter is measured 
along many directions. The absorption measurements are 
then inverted to yield the three dimensional distribution. 
Tomographic methods where first applied in astronomy by 
Bracewell d, for imaging microwave emitting regions on the 
surface of the Sun. The Radon transform forms the basis of 
the technique of Doppler tomography d which has been used 
to image accretion discs. 

We will now describe those aspects of the Radon transform, 
and its inversion, that are of immediate use to us. Readers can 
refer to an extensive collection of textbooks and monographs 
on the subject such as for an in-depth treatment of the 
subject. 


A. Definition 

Let x = (x\ ,X 2 , ■■ ■x n ) be a point in R". Let ^ and p define 
a hyperplane containing x, 

P = X ■ X = %\X! + %2*2 j \ - ZnXn- (1) 

The Radon transform, f(p, q ), of a function /(x) is defined as 

}\P-X) = / /(x) 8 (p - % ■ x)dx (2) 

Thus f(p , ) is the integral of /(x) over the hyperplane de¬ 
fined by p, c,. This integral is called a projection in the Radon 
transform literature. 

There exist several inversion methods for reconstructing 
/(x) from its radon transform f(p. ^ j. (In fact, an inversion 
method was provided by Radon in the same paper where the 
transform was first defined.) Here, we use the direct Fourier 
inversion method developed by Bracewell EE Cl which is 
based on a remarkable relation between the Fourier and Radon 
transforms of a function. 



FIG. 1: The relation between the Radon and Fourier transforms of a 
function. 

B. Relation to the Fourier Transform 

Let a point in Fourier space be denoted by, k = 
(k\. ko. ■ ■ ■ k N ). We adopt the following convention for the n- 
dimensional Fourier transform, &f, of a function fix), 

■m*) = J f(x)e- 2 ^dx (3) 

Let 4) be the one-dimensional Fourier transform, 

fj{p,1)e- 2 ™Pdp (4) 

It can be proved that, 

Jf(k) = ,^/OM). (5) 

That is, the n-dimensional Fourier transform of a function can 
be obtained by a suitable one-dimensional Fourier transform 
of its Radon transform. 

C. Direct Fourier Inversion Algorithm 

Eq . |5] is the basis of the direct Fourier inversion algorithm 
for reconstructing /(x) from its Radon transform f(p. c, ). The 
basic steps of this algorithm are as follows. 

1. Compute the one dimensional Fourier transform 

for a fixed k. 

2. Repeat the above step for all possible values of k to ob¬ 
tain the Fourier transform /(k) = ■'jj'ijk) in the Fourier 
domain. 

3. Compute the inverse Fourier transform, f to get 

/(x). 

Figure Q shows a schematic of the direct Fourier inversion 
method. 
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III. LISA RESPONSE TO THE BINARY POPULATION 
SIGNAL 


In this section, we set up much of our basic notation. We de¬ 
scribe, first, the response of LISA to a single binary and then 
use it to obtain the response to a binary population signal. We 
assume that a binary emits a monochromatic GW signal even 
though real binaries may not have circular orbits and, in addi¬ 
tion, will have some frequency evolution due to orbital decay 
by GW radiation. In some binaries, accretion from one star 
to another may drive frequency evolution. The assumption 
of a monochromatic signal should, however, be a fairly good 
approximation for most binaries at frequencies ~ 10 3 Hz or 
lower. In any case, it is the simplest model to use for a first 
study such as this. At low frequencies (< 5 mHz), where the 
binary population signal is expected, it is sufficient to use the 
LISA detector response given in | TUl . 


A. Signal from a single binary 



FIG. 2: The coordinate convention used to describe a source in the 
barycentric frame. 


We choose a Sun centered polar coordinate system with 
the orbital plane of the LISA centroid corresponding to po¬ 
lar angle 0 = 7i /2 and the azimuthal angle 0=0 corre¬ 
sponding to the position of LISA at the start, t = 0, of ob¬ 
servations. Let vv = ( sin 0 cos 0, sin 0 sin 0, cos0 ) be the 
unit vector pointing in the direction of source with angu¬ 
lar position (0,0). Define unit vectors 0 = dw/d6 and 
0 = (l/sin0)<9vv/<90 transverse to vv. The GW wave propa¬ 
gates in the —vv direction and the metric perturbation is given 

by. 


B. Detector Response 

The orbits of the three LISA satellites are designed such 
that the entire constellation will appear as a equilateral triangle 
in near rigid motion around the Sun. The plane of this triangle 
will be tilted at 60° to the orbital plane of its centroid and the 
triangle will rotate once in its plane during a year. We neglect 
corrections, such as the “flexing" of the whole configuration, 
to the basic picture above of LISA orbital motion. 

The signal phase observed at the LISA centroid will be 
doppler modulated due to orbital motion. 


fyuv(0 — h+(t) £+,/iv + h x (t) £x,jiv> 

e+ = 0 g) 0 — 0 (g) 0, 

£ x = 0(g)0 + 0(g)0. 


( 6 ) 

(7) 

( 8 ) 


<L(f) = f2f + DT> D (f)+4> o, (12) 

®D(t) = w-r(t) , (13) 

where r is the position vector of the centroid of LISA, 


where, h + (t) and h x ( t ) are the waveforms of the two indepen¬ 
dent polarizations of a GW in the Transverse-Traceless gauge. 


For a monochromatic binary source, the polarization wave¬ 
forms are given by. 


h+{t) 

h x (t) 

A+ 


A [A + cos2i//cos<t>(f) +A X sin2i//sin<t>(f)], (9) 
A [A x cos2t//sin<t>(f) — A + sin2y/cos<t>(f)] ,(10) 
1 + cos 2 e 


A x = cose 


( 11 ) 


where e is the angle between vv and the orbital angular mo¬ 
mentum vector of the binary, y/ denotes the angle between 
0 and the projection of the angular momentum vector on the 
plane transverse to vv, and A is an overall distance dependent 
factor. The signal phase observed in the Sun centered frame 
is <t>(f) = £lt + O 0 where £1 is twice the orbital frequency of 
the binary and 4>o is the initial phase at the beginning of the 
observation period. 


r = Rq (coso)gf,sinC00f,O) , (14) 

R . = 1 AU/c being the orbital radius of LISA, c the speed of 
light, and to. = 2 n/T @ , the orbital period corresponding to 
Tq = lyear. 

The detector response is given by the contraction of the sig¬ 
nal tensor, given by Eq. with the detector tensor, 

hu(t)=h flv :^ J v . (15) 

where the detector tensor 3>u is defined as follows. In the 
low frequency regime, where the armlength of the interferom¬ 
eter are small compared to the wavelength of the GW signal, 
the detector tensor can be written as the standard Michelson 
combination iTnll , 

9>u{t) = n I {t)®ni(t) — n J {t)®n J (t), (16) 

where nj represents the unit vectors along the 7 th arm of the 
interferometer. In the following we consider only one out of 
the two possible LISA interferometers. Hence, we drop the 
subscripts on hu(t) and S>u. 
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C. Response in terms of the Binary Distribution function 

Let p(0) be the number density of the binaries as func¬ 
tion of parameters 0 = ( 0 , 0 ,£ 2 , 17 ), where ?] is the subset 
(e, y/,A,<t>o). The number of binaries in the volume element 
d® at 0 is given by p (®)d®. We express p (0) as, 

p(0) = p a (£2|0,0,77)p'(0,0,77), (17) 

where p' is p marginalized over £2 and p_o is the conditional 
density of binaries in £2 for a given ( 0 , 0 , 17 ). 

For a given point 0, 0, I] in parameter space, the GW signal 
received at the origin of the Sun centered frame due to binaries 
in the volume element ddd(j)d 77 is 

h+,x(t',d,(j),ri) = f £/£2p a (£2|0,0,77)/z +jX (/;?7,£2) . (18) 


where Eq. [9] and M provide the expressions for /z + (/;T 7 ,£ 2 ) 
and /z x (/;/], £2) respectively. We can decompose 
h +]X (f; 0,0, 77 ) into a Fourier series. 


tv- 1 

A + , x (f m ;0,0,rj) = £ «''“*-*+,,<[*;0,0, 17 ] , (19) 

k =0 


where, for an observation period 7 G b s , the frequency a>k = 
2nk/T 0 \, s and t„, = zzzTobs/TV, N being the number of samples. 

The response, /*(/; 0,0, 77 ), of LISA to the GW signal com¬ 
ing from the volume element dQd(j)dri can be written as (see 
Eq.[j5}, 


/z(/; 0 , 0 , 77 ) = z/ 0 z/ 0 z/r 7 p / ( 0 , 0 , T 7 )/ 2 mv (/; 0 , 0 , 77 ) : , 


N -1 


V(f m ;0,0 ,j?) = £ ^ E e, ' fflA ' ((m+ ^ ( ' m ’ e ’ Ws «[^0’0^] h 

fl=+,x \i;=0 / 


r 


( 20 ) 

( 21 ) 


where we have included the doppler modulation <t>/) due to the 
orbital motion of LISA and £>(t) is given by Lq. 1 1 61 The total 
response of LISA to the binary population signal is. 




Ka(0,4>,t m ) 

K a [k;d,p} 

Ffl(0,0Tm) 


Jddd^dp h(t m \ 0 , 0 , rj) 
f dddtj) Y, T a ( 0 , 0 ,r,„) x 

J a=+,x 

7f a (0,0,f m + O o (t m ,0,0)) , (22) 

£^fc0,0]e io W", (23) 

J JT7p'(0,0,r])s a [A:;0,0,r7] , (24) 

e«,pv(0,0):^ v W (25) 


Now, 7y,(0,0,f) is a completely known function, as is 
<&£)(/, 0 , 0 ), and, hence, it can be replaced by a function 
G„(0,0,f + <!>,/(/, 0,0)) =F a (0,0,f). Thus, 


h{t m ) = JdOdip 7sr(0,0,f m + O D (h„, 0 , 0 )) , (26) 

where we have combined the integrand in Eq.|22]into a single 
function K(0. @.t). Since we do not know the binary distribu¬ 
tion function, K(6. @. t) is an unknown function. 


D. LISA Response as a Radon Transform 

It is clear from Eq. [25] that the response h(t) is obtained 
by integrating the function R( 0 , 0 ,/) on the surface ( 0 . 0 ./ + 
<f>£)(0,0,/)) parametrized by 0 and 0. Except for the fact that 


this surface is not a plane, h(t) is, therefore, nothing but the 
Radon transform (Sec. ED of the three dimensional function 
R(0,0,/). To cast h(t) as a standard Radon transform, we can 
make the following coordinate transformation which turns the 
surface of integration into a plane. Let x = (. x,y,z ), 


( y/2 R® sin 0 cos 0 
V2R @ sin 0 sin 0 

V2 (/ + <&/>) 

In terms of the new coordinates and using Eg. 1131 


(27) 


/ = —= — Rg sin 0 COS 0 COS 0 )g/ — 
v 2 

R® sin 0 sin 0 sin to 0 /, 


(28) 


and, hence, / = £ • x, where 

4 = —!p= ( — costdg/, —sin(d(.4. 1 ), (29) 

v 2 

In the new coordinates, 

h(t) = J K(x)8(t — 4 -x)rfx , (30) 


which casts the LISA response in the form of a standard 
Radon transform. (The Jacobian associated with the coordi¬ 
nate transformation has been absorbed in R(x).) 

Although the function K(x), Eq.[30| is not the same as the 
full distribution function p( 0 ), it already informs us about 
how the binary distribution appears in the three most inter¬ 
esting parameters, namely, 0, 0 and £2. In fact, for the sim¬ 
ple case of an isotropic detector antenna pattern the Discrete 
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Fourier Transform of K (x) along the z direction yields the dis¬ 
tribution of binaries in frequency Cl, for a given direction 8, 
(j), marginalized over the subset of parameters 77 . (This can be 
seen from Eq . 1231 and 1241 above.) In the following we will re¬ 
strict our attention to the reconstruction of K(x) and by taking 
its inverse Fourier transform along z, reconstruct the distri¬ 
bution of binaries along the spatial position (x,y) and signal 
frequency Cl. 


IV. TOMOGRAPHIC RECONSTRUCTION OF LISA 
BINARY DISTRIBUTION 


A (<C one year). The time series h(t) (the LISA response) 
can be divided into segments of length A to each of which a 
constant B, can be associated, such as B, (to) where to is the mid 
point of the segment. 

Under the above approximation, for a given B, (t) we now 
have more than one projection corresponding to (nearly) par¬ 
allel planes at distances t £ [to — A/2, to + A/2] from the ori¬ 
gin. (The sampling of the LISA response will give rise to a 
finite number of projections.) Hence, we can now use the di¬ 
rect Fourier inversion method. The steps involved in the latter 
were enumerated in Sec. mi Translated to the particular case 
here, they are as follows. 


Having shown that the response of LISA is a Radon trans¬ 
form of the binary distribution function (more precisely, its 
Fourier transform along the source frequency), we look at how 
and to what extent the latter can be reconstructed from the for¬ 
mer. As mentioned earlier there are several inversion methods 
for the Radon transform. In this paper we consider the direct 
Fourier inversion method that was outlined in Sec. im 

We use a simplified scenario for demonstrating reconstruc¬ 
tion: The response of LISA to a binary signal is assumed to 
involve only the orbital doppler frequency modulation and, 
although LISA will have three interferometers, the data is as¬ 
sumed to be the output of only one interferometer. This toy 
model is a conservative one since adding the amplitude and 
phase modulation due to the rotating antenna pattern of LISA 
and bringing in a second interferometer can only bring im¬ 
provements. 


A. Missing projections and approximate solution 

A Radon transform is invertible provided the integrals of 
the original function are known for all possible hyperplanes. 
The first problem we run into when inverting the LISA re¬ 
sponse is that we have, in fact, very few projections available 
to us. This is simply because of the fact that we have only one 
orbit producing the required doppler modulations for us. To 
construct more projections, we need multiple LISA detectors 
moving on different orbits around the Sun. 

The planes on which projections are available are specified 
by the unit normal 5 , defined in Eq. < 1291 . and the distance, t , of 
the plane from the origin. It is clear that / is confined to a two 
dimensional cone and does not cover all possible directions. 
Moreover, since ^ is a function oft, we have only one plane 
for any t whereas the Fourier inversion method (Sec HD re¬ 
quires that for a fixed ^projections be available for all possi¬ 
ble t. In the case of LISA, therefore, we have a severe shortage 
of projections and correspondingly our ability to reconstruct 
source position and frequency will be limited. 

It may appear that, given the above problem, the method of 
direct Fourier inversion is inapplicable to the inversion of the 
LISA response. However, we can take advantage of the fact 
that <^(f) changes very slowly, compared to the period of the 
sources LISA is designed to detect, to make an approximation 
that allows direct Fourier inversion to be used: Assume that 
B, (t ) is nearly constant over a sufficiently small time interval 


B. Numerical implementation of Direct Fourier Inversion 


First, the Discrete Fourier Transform of each segment 
(length A) of h(t ) is computed. Each segment is obtained by 
multiplying h(t) with a boxcar window function of length A. 
(We ignore end effects caused by the smaller length of the last 
few segments.) Consecutive segments have an overlap of half 
the segment length. The above step yields the the three di¬ 
mensional Fourier transform, K( k), of the function K(x) that 
is to be reconstructed (Fig.Q. 

The fact that B, is confined to a cone (Eq. 1291 restricts the 
Fourier wave-vector k to a cone too. Define polar coordinates 
(k. e, 1 //) in the Fourier space where e is the polar angle mea¬ 
sured from k z and ) f/ is the azimuthal angle measured from k x . 
Then k is confined to the cone defined by 


k 


k 

V2 


( - cos I//, 


— sin y/, 1 ) . 


(31) 


The next step is to compute the three dimensional inverse 
Fourier transform of K(k). Transforming to polar coordinates, 
the expression for the inverse Fourier transform can be written 
as, 


d\jjj ( 1 e-^ Kik ( xcos 'l'+ ysill v) x 

K(k,e = ^,\l/) ■ (32) 

Setting e = n/A accounts for the restriction of the first fourier 
transform to the cone defined in Eq. EU Having done the in¬ 
verse Fourier transform above, we will obtain K(x). 

Recall that in order to obtain the distribution, K(8,(j),Cl), 
of the binaries in frequency. Cl, the Fourier transform of K(x) 
along z needs to be computed for each (x.y) (which corre¬ 
spond to sky position). By definition, a doppler modulated 
monochromatic signal with carrier frequency fo originating 
at 8, 0 corresponds to a pure sinusoidal function of z with 
frequency fo/V2 at the corresponding location x, y. Let 
the Fourier transform along z of K(x) be g(k;x,y). Then 
G(k\x,y) g(k/V 2) is the desired source distribution along 
frequency, 

K(x) = J dke 2mkz g(k\x,y) J dke^ 2mkz G(t,x,y) (33) 


K(x) = jdk- 


a \fl Kikz 

7T 
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Thus, from Eq.|32]and|33] it follows that 

K(x,y,k) = (k, ip) . (34) 

is the distribution function in sky position and frequency. The 
integral over t// is approximated by a simple summation. (We 
have ignored overall constant factors in the above derivation.) 


C. Point Spread Function of tomographic reconstruction 

We have described the numerical implementation of the to¬ 
mographic reconstruction method above. The method is ex¬ 
plicitly linear, since it only involves two consecutive Fourier 
transforms. Consequently, it can be fully characterized in 
terms of its response to a delta function distribution, i.e., the 
point spread function (psf). 

The distribution function K(x) corresponding to a single 
monochromatic source is given by, 

^ S(x — a)8(y — b). (35) 

Here, (a, b) are the source location and T>o is the angular fre¬ 
quency of the signal. The Radon transform K(£,t) of the 
above source distribution is, 

K(lf) = JdxK{x)8(t-x-l)=e ia ° ( - , - a ^~ b ^ ) (36) 

where we have used the fact (c.f.. Sec. m that in the case 
of LISA, 4 is confined to a cone (c.f.. Eg. 12911. 

The first step in the inversion (c.f. Sec II V All is to break the 
time series h(t) = K(£(t),t) into sections of length A. The 
next step is to take the Fourier transform of each section. For 
the I th section, 

K(k{) = [ t,+A eim-ati-blh) e ~2mh dt , 

Jti 

= Ae~ ia °^ l+b ^ e 2m ’(/o-*)(<i+ A / 2 ) x 

sine ((/o — k)A) , (37) 

k, = k% , 


K(x) = exp 


(i&oz 

\V2 


that are integral multiples of 2tz/cQq, substituting K(k,e,yr) 
obtained above into Eq. 1341 gives 


K(x,y,k) 

pcosx 

psinx 

n 


2ni n Ak / TN /A x W 

— exp^,U-*)( I + —JJx 


sine (n (/o — k) A)J„ (2np) , 

(39) 

-2=(kx-f 0 a ), 

(40) 

2=(ky-f 0 b ), 

(41) 

2n(f 0 -k) 

5 

(42) 


where J„(x) is the Bessel function of the first kind of order 
»i. (There is no convenient representation for frequencies 
that are non-integral multiples of 2 7 t/&>q.) 

An important feature of the psf is the dependence (Eq. 
of the order n of J„ on frequency k. At the source frequency, 
the psf behaves like Jo(2nr) in the (x,y) plane, where r is 
the distance from the correct source location. However, when 
k /o, the behaviour is that of J n (x), n / 0. Now, the 
first maximum of J n (x) occurs at x ~ n. Therefore, at the 
correct source location, the psf falls to zero rapidly as one 
moves away from the source frequency. This fall off is much 
faster than the sine (n (/o — k) A) dependence. Thus, although 
naively one may expect a frequency resolution no better than 
~ 1/A due to the sectioning of the data, the frequency resolu¬ 
tion is actually ~ 1 /T a b s . 


V. EXPLORATORY RESULTS 


The psf described above has sidelobes that spread a point 
source distribution upon reconstruction. The sum of several 
sources, therefore, leads to maxima in the reconstructed den¬ 
sity which do not correspond to a real source. It should be 
possible to obtain a cleaner reconstruction by applying sophis¬ 
ticated deconvolution methods such as CLEAN 01. However, 
this is a fairly non-trivial task which needs to be addressed 
carefully. We leave this to a future work and report only ex¬ 
ploratory results in this paper. 


A. LISA response and reconstruction parameters 


where | (/,- + A/2). Finally, the three dimensional inverse 
Fourier transform is computed as described in Sec. II V Bl In 
polar coordinates, 

K(k,e=^,\l?j = A e-i^Macosy+bsiny) x 

ex p( 2 »(/o —*)(Y; + §))x 
sine ((/o ~k)A) , (38) 

where we have used the fact that yt = co t, ft) Q being the or¬ 
bital period of LISA (c.f., Eq. 1291) . For source frequencies 


The time series h(t) representing the response of LISA to 
the binary population signal is generated as follows. Doppler 
modulated sinusoidal signals, 

s j (t)=Acos(Q.j(t-Xj£ l (t)-y j Z 2 (t))) , (43) 

are generated for a set of randomly chosen sky locations 
(x,y) (c.f., Eq. |^7] ) and frequencies £2. The spacing A k be¬ 
tween consecutive values of k = Q./2tz values is chosen to be 
bmm$k < Ak < b mm 8k, where 8k = l/T^bs, 7ob s being the ob¬ 
servation period (taken to be 1 year). (The examples shown 
here correspond to different choices for the integers b m j n and 
^rnax-) We ignore the initial phase of the signal because we 
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do not pursue reconstruction of the distribution in parameters 
beyond (x,y,Cl). The binary population signal is, 

h( f) = E‘ s f( f ) ■ (44) 

j 

We take / = 1..... 100 and k\ = 2 mHz. The distribution func¬ 
tion is then reconstructed only in a range 32 x 8k centered at 
a frequency > 10 x 8k + k\. This procedure accounts for the 
influence of binaries lying outside the frequency band of im¬ 
mediate interest. All binary signals have the same matched 
filtering signal to noise ratio (snr) p, defined for white noise 
as. 



For the reconstruction, we split h(t) into segments that are A = 
1 day long (c. f. SecEV} with an overlap between consecutive 
segments of half this length. The sampling frequency of the 
data is set at 0.1 Hz. 


B. Reconstruction examples 

In the following, we use \K(x.y,k)\ (Ea. l34> to study tomo¬ 
graphic reconstruction. We use a simple method for identify¬ 
ing sources: First, \K(x,y,k)\ is obtained for LISA response 
consisting of only white, Gaussian noise. From this, we es¬ 
timate a value t] such that it is the smallest supremum over 
all values of \K(x,y,k)\ in a specified range. Each sky map 
\K(x,y,k)\, obtained for data containing both noise and the bi¬ 
nary population signal h(t), is then normalized by 77 . (The 
same noise realization is used in both the noise only and noise 
plus signal reconstructions.) A threshold is then applied to the 
normalized sky maps and any isolated cluster of pixels in the 
(x,y) plane with values greater than the threshold is taken to 
be an individual source. 

a. Noise only - Fig. 0 shows |/f(x,y,k)| at two widely 
separated frequencies k\ and ki, ki — k\ = 13 x 8k, when h(t) 
is a realization of white, Gaussian noise. From the values of 
\K(x,y,kj)\ over the bandwidth kj £ [k\ + 10 x 8k,k\ +42 x 
8k] mentioned above, we estimate 77 , the smallest supremum 
of all pixel values in this frequency band (and all values of 
(x,y)) to be 77 = 3 x 10 4 . (The overall scale is arbitrary here.) 
The effect of the psf is apparent from the presence of slow 
variations in the amplitude of the reconstructed distribution. 

b. Noise plus signal from well isolated, strong binaries - 

This example demonstrates the recovery of fairly large snr 

binaries that are well isolated in frequency. Fig. 0 and 0 show 
\K{x,y,k)\/r) for a realization of LISA response that has a bi¬ 
nary population signal h(t) added to the same noise realiza¬ 
tion as Fig.0 For the signal we choose p = 15, b mm = 10 and 
^max = 15. 

After normalization with 77 , a threshold of \K[x,y,k)\/ r\ = 
4 is applied to all the sky maps. The contours of \K(x,y,k)\/f] 
greater than this threshold are superimposed on the sky maps. 
We find that only the maximia associate with the true source 
locations cross the threshold. (There are two minor peaks that 


cross the threshold but are not associated with the main lobe 
of the psf.) The maxima of \K(x,y,k)\ do not coincide exactly 
with the true source location. This is due to the fact that the 
frequencies of the sources are not exactly coincident with the 
those associated with the maps. (The binary frequencies need 
not correspond exactly to integral multiples of 8k. When this 
is the case, we show the source location in the two sky maps 
on either side of the true source frequency.) 

Another feature to note in the figures is the rapid falloff of 
the psf at the source location with a shift of one frequency bin 
from the source frequency. This behavior is expected from the 
analytical form given in Sec. HV Cl 

c. Noise plus signal from well isolated, weak binaries - 

Fig- 0 and |7| show \K{x,y,k)\/r\ for a realization of LISA 

response that has a binary population signal h(t) added to the 
same noise realization as Fig. 0 For the signal we choose 
P = 7, b m i n = 10 and Z? max = 15. 

After normalization with 77 , a threshold of \K{x.,y,k)\/r\ = 
1.8 is applied to all the sky maps. The contours of 
\K(x,y,k)\/r) greater than this threshold are superimposed on 
the sky maps. We find that, for this threshold, the maxima as¬ 
sociated with the true source locations cross the threshold. A 
few false peaks emerge due to the noise. This example pro¬ 
vides an estimate of the minimum signal strength required for 
a probability of detection near unity at a rate of false source 
detection that is comparable to the rate of occurrence of true 
binaries. 

d. Noise plus signal from close, weak binaries - Fig.[S] 
and|9]show \K(x,y,k)\/i] for a realization of LISA response 
that has a binary population signal h(t) added to the same 
noise realization as Fig. 0 For the signal we choose p = 7, 
bmin = 4 and b max = 8 . After normalization with 77 , a thresh¬ 
old of \K(x,y,k)\/r) = 1.8 is applied to all the sky maps. The 
contours of \K(x,y,k)\/rj greater than this threshold are super¬ 
imposed on the sky maps. 

In this example we see the strong effect of the sidelobes of 
the psf. The number of false sources increases compared to 
the earlier example. However, the true source locations con¬ 
tinue to have strong maxima above the threshold. A deconvo¬ 
lution method should be used in this situation to mitigate the 
effect of the sidelobes. 


VI. CONCLUSIONS 

We present the first steps of a tomographic approach to the 
reconstruction of the LISA Galactic binary distribution. The 
connection between the LISA response and the Radon trans¬ 
form of the binary distribution function is established. The 
direct Fourier inversion method is then applied to invert the 
response and reconstruct the distribution function. 

A full characterization of the performance of the method 
will require Monte Carlo simulations with realistic LISA data. 
In this paper, we have presented exploratory results, obtained 
using a simple source identification procedure. We find that 
near 2 mHz, it is possible to recover binary sources having an 
snr of about 7 and separated by about 10 (one year) frequency 
bins. This is close to the expected density of binary sources 
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FIG. 4: Sky maps of \K(x,y,k)\/r\ for T) = 3 x 1(L 4 , p = 15 and sources set apart by a minimum of 10 frequency bins. One frequency 
bin has size 8f = 3.17 x 10~ 8 corresponding to 1 /7obs where T 0 \, s = 1 year is the period of observation. Contours are drawn at the levels 
[4,4.2,4.4,4.6,4.8, 5.0]. The frequency kj corresponding to the panel in row 1 and column m, from top and left respectively, is given by 
Aq,i + ((/ — 1) x 4 + w) x 8f. The location of a source is indicated in the title of the panels that are closest in frequency on either side of the 
source frequency. The symbol marks the source location. 


around this frequency. For a higher density of sources, one 
must contend with a large number of false sources that arise 
from the superposition of the sidelobes of the psf. Deconvo¬ 
lution methods such as CLEAN should be useful in reducing 
this effect. 

The reconstruction method presented in the paper does not 


take into account (i) the amplitude and phase modulations due 
to the rotating antenna pattern of LISA, and (ii) the presence 
of two interferometers in the LISA configuration. Inclusion of 
these factors into the method is straightforward and doing so 
should lead to better performance. 

Since the nature of the orbit of LISA does not furnish a 
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FIG. 5: Continuation of Fig. [4] 


sufficient number of projections for exact inversion, we had 
to take recourse to a simple approximation in order to apply 
the Direct Fourier inversion method. This situation of missing 
projections is encountered quite often in real world imaging 
applications, though probably not at the severe level of LISA. 
Thus, a better approach to inversion may be possible than the 
straightforward one used here. Work on this and the issues 
mentioned above is in progress. 

The attractive features of the tomographic approach are the 
low and fixed computational cost of the Radon inverse trans¬ 
form and the linearity of the method. For 1 year of data 
sampled at 0.1 Hz, a typical Pentium-4 desktop needs about 
2 hours to produce the reconstructions shown in this paper. 
This is independent of the number of sources but does de¬ 
pend on the number of frequencies at which a sky map is ob¬ 


tained. However, the method can be parallelized quite easily 
and, hence, computational cost does not appear to be a signif¬ 
icant issue. 
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FIG. 6: Sky maps of \K{x,y,K)\/r] for r\ = 3 x 10~ 4 , p = 7 and sources set apart by a minimum of 10 frequency bins. One frequency bin has 
size Sf = 3.17 x 10~ 8 corresponding to 1/T 0 \, S where r o b s = 1 year is the period of observation. Contours are drawn at the levels [1.8,1.9,2.0], 
The frequency kj corresponding to the panel in row / and column m, from top and left respectively, is given by k\^ + ((/ — 1) x 4 + /;i) x Sf. 
The location of a source is indicated in the title of the panels that are closest in frequency on either side of the source frequency. The 
symbol marks the source location. 
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FIG. 7: Continuation of Fig.[6]in frequency. 
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FIG. 8: Sky maps of \K(x,y,k)\/r\ for T) = 3 x 10~ 4 , p = 7 and sources set apart by a minimum of 4 frequency bins. One frequency bin has 
size 8f = 3.17 x 10~ 8 corresponding to l/r o b s where T 0 \, s = 1 year is the period of observation. Contours are drawn at the levels [1.8,1.9,2.0]. 
The frequency kj corresponding to the panel in row / and column m, from top and left respectively, is given by Aq i + ((/ — 1) x 4 + ;;;) x 8f. 
The location of a source is indicated in the title of the panels that are closest in frequency on either side of the source frequency. The 
symbol marks the source location. 
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FIG. 9: Continuation of Fig.[8] 



































































